# Импортируем необходимые модули
import numpy as np
import matplotlib.pyplot as plt
import plotly
import plotly.graph_objs as go
import plotly.express as px
from plotly.subplots import make_subplots
from ipywidgets import interact
import warnings
warnings.filterwarnings("ignore")
# Определяем функцию, вычисляющую правую часть системы
def sir_model(t, x, N, beta, gamma):
S, I, R = x
dS_dt = -beta*I*S/N
dI_dt = beta*I*S/N - gamma*I
dR_dt = gamma*I
xdot = np.array([dS_dt, dI_dt, dR_dt])
return xdot
# Задаём параметры и начальные условия
beta = 0.32
gamma = 0.12
N = 200000
I0 = 100
R0 = 0
S0 = N - I0 - R0
t0 = 0
tf = 100
dt = 0.01
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
# Реализуем метод Рунге-Кутты 4-го порядка
"""
t0 и tf - начальный и конечный момент времени
dt - шаг времени
x0 - массив значений в начальный момент времени
"""
def RungeKutta4(f, x0, t0, tf, dt):
# Массив всех t от t0 до tf с шагом dt
t = np.arange(t0, tf, dt)
nt = t.size
# Массив из nx строк и nt столбцов, пока что заполненный нулями
nx = x0.size
x = np.zeros((nx, nt))
# Заполняем первый столбец (при t = 0) начальными значениями
x[:,0] = x0
for k in range(nt - 1):
# Коэфициенты для метода РК
k1 = f(t[k], x[:,k])
k2 = f(t[k] + dt/2, x[:,k] + dt*k1/2)
k3 = f(t[k] + dt/2, x[:,k] + dt*k2/2)
k4 = f(t[k] + dt, x[:,k] + dt*k3)
dx = dt*(k1 +2*k2 +2*k3 +k4)/6
# Запись новых значений в текущий момент времени
x[:, k+1] = x[:,k] + dx
return x, t
# Решение методом Рунге-Кутты 4-го порядка
x, t = RungeKutta4(f, x0, t0, tf, dt)
# Визуализируем решения для метода Рунге-Кутты 4-го порядка
def plot(x, t):
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=x[0,:], name="Восприимчивые (S)"))
fig.add_trace(go.Scatter(x=t, y=x[1,:], name="Заражённые (I)"))
fig.add_trace(go.Scatter(x=t, y=x[2,:], name="Выздоровевшие (R)"))
fig.update_layout(legend_orientation="h", xaxis_title="Время (t)", margin=dict(l=0, r=0, t=30, b=0))
fig.show()
fig1 = go.Figure(data = [go.Scatter3d(x=x[0,:], y=x[1,:], z=x[2,:])])
fig1.update_layout(scene = dict(
xaxis_title="Восприимчивые (S)",
yaxis_title="Заражённые (I)",
zaxis_title="Выздоровевшие (R)"),
width=700,
margin=dict(l=0, r=0, t=30, b=0))
fig1.show()
plot(x, t)
# Добавим интерактив для метода Рунге-Кутты 4-го порядка
def interact_rk4(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100):
dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x, t = RungeKutta4(f, x0, t0, tf, dt)
plot(x, t)
interact(interact_rk4, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(0, N, 10), tf=(10, 1000, 1))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_rk4(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100)>
# Реализуем метод Рунге-Кутты 4-го порядка с правилом 3/8
"""
t0 и tf - начальный и конечный момент времени
dt - шаг времени
x0 - массив значений в начальный момент времени
"""
def RungeKutta_3_8(f, x0, t0, tf, dt):
# Массив всех t от t0 до tf с шагом dt
t = np.arange(t0, tf, dt)
nt = t.size
# Массив из nx строк и nt столбцов, пока что заполненный нулями
nx = x0.size
x = np.zeros((nx, nt))
# Заполняем первый столбец (при t = 0) начальными значениями
x[:,0] = x0
for k in range(nt - 1):
# Коэфициенты для метода РК
k1 = f(t[k], x[:,k])
k2 = f(t[k] + dt/3, x[:,k] + dt*k1/3)
k3 = f(t[k] + 2*dt/3, x[:,k] - dt*k1/3 + dt*k2)
k4 = f(t[k] + dt, x[:,k] + dt*k1 - dt*k2 +dt*k3)
dx = dt*(k1 +3*k2 +3*k3 +k4)/8
# Запись новых значений в текущий момент времени
x[:, k+1] = x[:,k] + dx
return x, t
# Решение методом Рунге-Кутты 4-го порядка с правилом 3/8
x, t = RungeKutta_3_8(f, x0, t0, tf, dt)
# Визуализируем решение методом Рунге-Кутты 4-го порядка с правилом 3/8
plot(x, t)
# Добавим интерактив для метода Рунге-Кутты 4-го порядка с правилом 3/8
def interact_rk38(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100):
dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x, t = RungeKutta_3_8(f, x0, t0, tf, dt)
plot(x, t)
interact(interact_rk38, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_rk38(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100)>
# Реализуем метод Ральстона 3-го порядка
"""
t0 и tf - начальный и конечный момент времени
dt - шаг времени
x0 - массив значений в начальный момент времени
"""
def Rolstone3(f, x0, t0, tf, dt):
# Массив всех t от t0 до tf с шагом dt
t = np.arange(t0, tf, dt)
nt = t.size
# Массив из nx строк и nt столбцов, пока что заполненный нулями
nx = x0.size
x = np.zeros((nx, nt))
# Заполняем первый столбец (при t = 0) начальными значениями
x[:,0] = x0
for k in range(nt - 1):
# Коэфициенты для метода
k1 = f(t[k], x[:,k])
k2 = f(t[k] + dt/2, x[:,k] + dt*k1/2)
k3 = f(t[k] + 3*dt/4, x[:,k] + 3*dt*k2/4)
dx = dt*(2*k1/3 + k2 + 4*k3/3)/3
# Запись новых значений в текущий момент времени
x[:, k+1] = x[:,k] + dx
return x, t
# Решение методом Ральстона 3-го порядка
x, t = Rolstone3(f, x0, t0, tf, dt)
# Визуализируем решение методом Ральстона 3-го порядка
plot(x, t)
# Добавим интерактив для метода Ральстона 3-го порядка
def interact_r3(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100):
dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x, t = Rolstone3(f, x0, t0, tf, dt)
plot(x, t)
interact(interact_r3, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_r3(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100)>
"""
t0 и tf - начальный и конечный момент времени
dt - шаг времени
x0 - массив значений в начальный момент времени
"""
def Hoyne3(f, x0, t0, tf, dt):
# Массив всех t от t0 до tf с шагом dt
t = np.arange(t0, tf, dt)
nt = t.size
# Массив из nx строк и nt столбцов, пока что заполненный нулями
nx = x0.size
x = np.zeros((nx, nt))
# Заполняем первый столбец (при t = 0) начальными значениями
x[:,0] = x0
for k in range(nt - 1):
# Коэфициенты для метода
k1 = f(t[k], x[:,k])
k2 = f(t[k] + dt/3, x[:,k] + dt*k1/3)
k3 = f(t[k] + 2*dt/3, x[:,k] + 2*dt*k2/3)
dx = dt*(1*k1/4 + 0*k2 + 3*k3/4)
# Запись новых значений в текущий момент времени
x[:, k+1] = x[:,k] + dx
return x, t
x, t = Hoyne3(f, x0, t0, tf, dt)
plot(x, t)
def interact_h3(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100):
dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x, t = Hoyne3(f, x0, t0, tf, dt)
plot(x, t)
interact(interact_h3, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_h3(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100)>
def ssprk3(f, x0, t0, tf, dt):
# Массив всех t от t0 до tf с шагом dt
t = np.arange(t0, tf, dt)
nt = t.size
# Массив из nx строк и nt столбцов, пока что заполненный нулями
nx = x0.size
x = np.zeros((nx, nt))
# Заполняем первый столбец (при t = 0) начальными значениями
x[:,0] = x0
for k in range(nt - 1):
# Коэфициенты для метода
k1 = f(t[k], x[:,k])
k2 = f(t[k] + dt, x[:,k] + dt*k1)
k3 = f(t[k] + dt/2, x[:,k] + dt*k1/4 + dt*k2/4)
dx = dt*(k1/2 + k2/2 + 2*k3)/3
# Запись новых значений в текущий момент времени
x[:, k+1] = x[:,k] + dx
return x, t
x, t = ssprk3(f, x0, t0, tf, dt)
plot(x, t)
def interact_ssprk3(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100):
dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x, t = ssprk3(f, x0, t0, tf, dt)
plot(x, t)
interact(interact_ssprk3, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_ssprk3(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100)>
def compare_solutions(t0, tf, dt, x1, x2):
nx = x0.size
t = np.arange(t0, tf, dt)
nt = t.size
x = np.zeros((nx, nt))
for k in range(nt - 1):
x[:, k] = x1[:,k] - x2[:,k]
return x, t
def get_dt(func, acc = 0.00001, beta=0.32, gamma=0.12, N=200000, I0=100, tf=100):
R0 = 0
S0 = N - I0 - R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
dt = tf/2-1
flag1 = 0
flag2 = 0
while dt >= 0.00001:
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = func(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
for k in range(t.size - 1):
for m in range(len(x[:, k])):
if abs(x[:, k][m]) / x1[:, k][m] > acc:
flag1 = 1
break
if flag1:
flag1 = 0
flag2 = 1
break
if not flag2:
print("Необходимый максимальный шаг:", dt, "ед.")
return dt
else:
flag2 = 0
dt -= 0.01
def get_kl(func, beta=0.32, gamma=0.12, N=200000, I0=100, tf=100):
acc = float(input("Введите необходимую точность: "))
dt = get_dt(func, acc)
t = np.arange(t0, tf, dt)
nt = t.size
print("Необходимо рассчитать правую часть системы раз:",nt)
def interact_comp1(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100, dt = 0.01):
# dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = RungeKutta_3_8(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
plot(x, t)
interact(interact_comp1, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1), dt=(0.01, 1, 0.01))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_comp1(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100, dt=0.01)>
def interact_comp1_text(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100, dt = 0.01):
# dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = RungeKutta_3_8(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
for k in range(t.size - 1):
print(x[:, k], t[k])
interact(interact_comp1_text, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1), dt=(0.01, 1, 0.01))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_comp1_text(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100, dt=0.01)>
get_kl(RungeKutta_3_8)
Введите необходимую точность: 0.0001 Необходимый максимальный шаг: 2.870000000001161 ед. Необходимо рассчитать правую часть системы раз: 35
plot(x, t)
def interact_comp2(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100, dt = 0.01):
# dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = ssprk3(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
plot(x, t)
interact(interact_comp2, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1), dt=(0.01, 1, 0.01))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_comp2(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100, dt=0.01)>
def interact_comp2_text(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100, dt = 0.01):
# dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = ssprk3(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
for k in range(t.size - 1):
print(x[:, k], t[k])
interact(interact_comp2_text, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1), dt=(0.01, 1, 0.01))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_comp2_text(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100, dt=0.01)>
get_kl(ssprk3)
Введите необходимую точность: 0.0001 Необходимый максимальный шаг: 0.4300000000011781 ед. Необходимо рассчитать правую часть системы раз: 233
def interact_comp3(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100, dt = 0.01):
# dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = Hoyne3(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
plot(x, t)
interact(interact_comp3, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1), dt=(0.01, 1, 0.01))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_comp3(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100, dt=0.01)>
def interact_comp3_text(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100, dt = 0.01):
# dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = Hoyne3(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
for k in range(t.size - 1):
print(x[:, k], t[k])
interact(interact_comp3_text, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1), dt=(0.01, 1, 0.01))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_comp3_text(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100, dt=0.01)>
get_kl(Hoyne3)
Введите необходимую точность: 0.0001 Необходимый максимальный шаг: 0.4400000000011781 ед. Необходимо рассчитать правую часть системы раз: 228
def interact_comp4(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100, dt = 0.01):
# dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = Rolstone3(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
plot(x, t)
interact(interact_comp4, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1), dt=(0.01, 1, 0.01))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_comp4(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100, dt=0.01)>
def interact_comp4_text(beta=0.32, gamma=0.12, N=200000, I0=100, tf = 100, dt = 0.01):
# dt = 0.01
R0=0
S0 = N-I0-R0
t0 = 0
x0 = np.array([S0, I0, R0])
f = lambda t, x: sir_model(t, x, N, beta, gamma)
x1, t1 = RungeKutta4(f, x0, t0, tf, dt)
x2, t2 = Rolstone3(f, x0, t0, tf, dt)
x, t = compare_solutions(t0, tf, dt, x1, x2)
for k in range(t.size - 1):
print(x[:, k], t[k])
interact(interact_comp4_text, beta=(0, 2, 0.02), gamma=(0, 0.5, 0.01), N=(100, 400000, 100), I0=(10, N, 10), tf=(10, 1000, 1), dt=(0.01, 1, 0.01))
interactive(children=(FloatSlider(value=0.32, description='beta', max=2.0, step=0.02), FloatSlider(value=0.12,…
<function __main__.interact_comp4_text(beta=0.32, gamma=0.12, N=200000, I0=100, tf=100, dt=0.01)>
get_kl(Rolstone3)
Введите необходимую точность: 0.0001 Необходимый максимальный шаг: 0.4300000000011781 ед. Необходимо рассчитать правую часть системы раз: 233